!p.multi=[0,2,1]
!p.charsize=1.5
restore,'solardata.sav'
l=64
rr=6.96e8      ;m
rmin=0.6
rmax=0.96
rds=rmin*rr     ;m
rdsi=1./rds
GG=6.674e-11   ;N(m/kg)^2
M=1.988e30     ;kg
dz=(rmax-rmin)*rr/double(l-1)
r=rds+indgen(l)*dz
g=GG*M/r^2     ;m/s^2
g0=g[0]
g0=710.
z = indgen(l-1)*dz

; bottom values

;rho00=  420.920
;p00=   1.71889e+13
;T00=   2.97382e+06
;th00= 1.14*T00       ;2.06008e+06
rho00= 506.84
p00=   2.1404365e+14
T00=   3.12697e6
th00= 1.1*T00       ;2.06008e+06

rg=13732.
cp=2.5*rg
cap=rg/cp
capi=1./cap
ss=g0/(cp*th00)

rho = rho00*(1.-(ss*z)/(1.+z*rdsi))^(capi-1)

;pm0=p00*(1.-cap*z*sdi)**capi ; CAUTION: max. alt. exists
;tm0=T00*(1.-cap*z*sdi)       ; CAUTION: max. alt. exists

plot,r/rr,rho
oplot,re,rho_new,li=2

st=-3./(th00*(l-1)*dz)
;the=th00*(1.+ st*z)
t0=t00*(1.-ss*z/(1.+z*rdsi))

plot,r/rr,t0,yr=[min(t0),max(t0)]
oplot,re,t_new,li=2

!p.multi=0
end

